/* ***** BEGIN LICENSE BLOCK *****
   Version: Apache 2.0/GPL 3.0/LGPL 3.0

   CCT - Computational Chemistry Tools
   JMolEditor - Java Molecules Editor

   Copyright 2008 Dr. Vladislav Vasilyev

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.

   Contributor(s):
     Dr. Vladislav Vasilyev <vvv900@gmail.com>       (original author)

  Alternatively, the contents of this file may be used under the terms of
  either the GNU General Public License Version 2 or later (the "GPL"), or
  the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  in which case the provisions of the GPL or the LGPL are applicable instead
  of those above. If you wish to allow use of your version of this file only
  under the terms of either the GPL or the LGPL, and not to allow others to
  use your version of this file under the terms of the Apache 2.0, indicate your
  decision by deleting the provisions above and replace them with the notice
  and other provisions required by the GPL or the LGPL. If you do not delete
  the provisions above, a recipient may use your version of this file under
  the terms of any one of the Apache 2.0, the GPL or the LGPL.

 ***** END LICENSE BLOCK *****/

package cct.gaussian;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.RandomAccessFile;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.List;
import java.util.Set;
import java.util.StringTokenizer;
import java.util.logging.Logger;

import cct.Constants;
import cct.interfaces.AtomInterface;
import cct.interfaces.MoleculeInterface;
import cct.modelling.ChemicalElements;
import cct.modelling.MolecularGeometry;
import cct.modelling.StructureManagerInterface;
import cct.vecmath.vPoint3f;

/**
 * <p>Title: Computational Chemistry Tookit</p>
 *
 * <p>Description: </p>
 *
 * <p>Copyright: Copyright (c) 2005-2010 Dr. Vladislav Vassiliev</p>
 *
 * <p>Company: ANU</p>
 *
 * @author Dr. Vladislav Vassiliev
 * @version 1.0
 */
public class G03MDTrajectory
    implements StructureManagerInterface {

   public static final String TOTAL_ENERGY_KEY = "Total energy (au)";
   private static final String STANDARD_GEOMETRY_KEY = "MD";

   private boolean debug = true;
   private List<GaussianAtom> gAtoms = null;
   private List<GaussianSnapshot> snapshots = new ArrayList<GaussianSnapshot> ();
   private Set<String> propertiesToChart = new HashSet<String> ();
   static final Logger logger = Logger.getLogger(G03MDTrajectory.class.getCanonicalName());

   public G03MDTrajectory() {
   }

   /**
    * Contrary to any similar subroutines it returns the FIRST snapshot
    * @param molec MoleculeInterface -  molecule
    * @throws Exception
    */
   public void getMolecularInterface(MoleculeInterface molec) throws
       Exception {
      if (molec == null) {
         throw new Exception(this.getClass().getCanonicalName() + " : molec == null");
      }
      if (gAtoms.size() < 1) {
         throw new Exception(this.getClass().getCanonicalName() + " : no Gaussian atoms");
      }

      molec.addMonomer("GAUSS");

      for (int i = 0; i < gAtoms.size(); i++) {
         GaussianAtom ga = gAtoms.get(i);
         AtomInterface atom = molec.getNewAtomInstance();
         if (ga.name == null || ga.name.length() < 1) {
            atom.setName(ChemicalElements.getElementSymbol(ga.getElement()));
         }
         else {
            atom.setName(ga.name);
         }
         atom.setAtomicNumber(ga.getElement());
         atom.setXYZ(ga.xyz[0], ga.xyz[1], ga.xyz[2]);

         molec.addAtom(atom);
      }

   }

   public String[] getAvailPropToChart() {
      if (propertiesToChart.size() < 1) {
         return null;
      }
      String[] sa = new String[propertiesToChart.size()];
      propertiesToChart.toArray(sa);
      return sa;
   }

   public MoleculeInterface parseMolecularGeometryOnly(MoleculeInterface mol, String filename) throws Exception {
      try {
         BufferedReader in = new BufferedReader(new FileReader(filename));
         return parseMolecularGeometryOnly(in, mol);
      }
      catch (Exception ex) {
         throw new Exception("parseMolecularGeometryOnly: ERROR opening file " + filename);
      }
   }

   public MoleculeInterface parseMolecularGeometryOnly(BufferedReader in, MoleculeInterface mol) throws Exception {
      /*
                                  Trajectory Starting Point:
        ---------------------------------------------------------------------
        Center     Atomic     Atomic              Coordinates (Angstroms)
        Number     Number      Type              X           Y           Z
        ---------------------------------------------------------------------
           1          6      20001001       -0.241185   -0.271881   -0.382124
           2          1      20001035       -1.177106   -0.664240   -0.014125
           3          1      20001035       -0.315056    0.809277   -0.444002
           4          1      20001035       -0.048133   -0.668818   -1.370969
           5          8      20001022        0.768355   -0.655283    0.526139
           6          1      20001031        1.633769   -0.371391    0.194276
           7          8      20001021       -4.053620   -1.712267   -3.472290
           8          1      20001030       -4.735666   -1.205766   -3.951382
           9          1      20001030       -4.373260   -2.632389   -3.570404
          10          8      20001021       -1.639883   -2.524523   -2.527601
          11          1      20001030       -0.987622   -1.996276   -3.026744
          12          1      20001030       -2.490365   -2.190708   -2.882675
       ....
         424          8      20001021        4.214621    4.206208    5.313912
         425          1      20001030        3.628498    4.920676    5.625427
         426          1      20001030        4.840445    4.696562    4.738038
        ---------------------------------------------------------------------

       */
      String line;
      try {

         // --- Finding "Trajectory Starting Point"
         while ( (line = in.readLine()) != null) {
            if (line.contains("Trajectory Starting Point:")) {
               break;
            }
         }

         if (line == null) {
            throw new Exception("parseMolecularGeometryOnly: ERROR: Unxpected End of file while looking for the start of data");
         }

         // --- Skip 4 lines
         for (int i = 0; i < 4; i++) {
            line = in.readLine();
            if (line == null) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Unxpected End of file while looking for the start of coordinates");
            }
         }

         // --- Reading atoms
         //    1          6      20001001       -0.241185   -0.271881   -0.382124

         while ( (line = in.readLine()) != null) {
            if (line.contains("-------------------------")) {
               break;
            }

            StringTokenizer st = new StringTokenizer(line, " ");
            if (st.countTokens() != 6) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: coordinate input has 6 tokens\n  Got: " + line);
            }

            AtomInterface atom = mol.getNewAtomInstance();

            // --- Getting element

            // --- Skip 1 token

            st.nextToken();

            int element;
            try {
               element = Integer.parseInt(st.nextToken());
               if (element < 1) {
                  element = 0;
               }
               atom.setAtomicNumber(element);
               atom.setName(ChemicalElements.getElementSymbol(element));
            }
            catch (Exception ex) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Error while parsing atom: " +
                   (mol.getNumberOfAtoms() + 1) + " Cannot parse atom number: " + line);
            }

            // --- Getting atomic type

            // --- Skip it for now

            st.nextToken();

            // --- Getting x,y,z

            float xyz;
            try {
               xyz = Float.parseFloat(st.nextToken());
               atom.setX(xyz);
               xyz = Float.parseFloat(st.nextToken());
               atom.setY(xyz);
               xyz = Float.parseFloat(st.nextToken());
               atom.setZ(xyz);
            }
            catch (Exception ex) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Error while parsing atom: " +
                   (mol.getNumberOfAtoms() + 1) + " Cannot parse atom's coordinate(s): " + line);
            }

            mol.addAtom(atom);

         }

         in.close();
      }
      catch (Exception ex) {
         throw new Exception("parseMolecularGeometryOnly: ERROR: " + ex.getMessage());
      }

      return mol;
   }

   static public List<GaussianAtom> parseMolecularGeometryOnly(RandomAccessFile raf) throws Exception {
      List<GaussianAtom> atoms = new ArrayList<GaussianAtom> ();
      /*
                                Trajectory Starting Point:
              ---------------------------------------------------------------------
              Center     Atomic     Atomic              Coordinates (Angstroms)
              Number     Number      Type              X           Y           Z
              ---------------------------------------------------------------------
         1          6      20001001       -0.241185   -0.271881   -0.382124
         2          1      20001035       -1.177106   -0.664240   -0.014125
         3          1      20001035       -0.315056    0.809277   -0.444002
         4          1      20001035       -0.048133   -0.668818   -1.370969
         5          8      20001022        0.768355   -0.655283    0.526139
         6          1      20001031        1.633769   -0.371391    0.194276
         7          8      20001021       -4.053620   -1.712267   -3.472290
         8          1      20001030       -4.735666   -1.205766   -3.951382
         9          1      20001030       -4.373260   -2.632389   -3.570404
        10          8      20001021       -1.639883   -2.524523   -2.527601
        11          1      20001030       -0.987622   -1.996276   -3.026744
        12          1      20001030       -2.490365   -2.190708   -2.882675
             ....
       424          8      20001021        4.214621    4.206208    5.313912
       425          1      20001030        3.628498    4.920676    5.625427
       426          1      20001030        4.840445    4.696562    4.738038
              ---------------------------------------------------------------------

       */
      String line;
      try {

         // --- Finding "Trajectory Starting Point"
         while ( (line = raf.readLine()) != null) {
            if (line.contains("Trajectory Starting Point:")) {
               break;
            }
         }

         if (line == null) {
            throw new Exception("parseMolecularGeometryOnly: ERROR: Unxpected End of file while looking for the start of data");
         }

         // --- Skip 4 lines
         for (int i = 0; i < 4; i++) {
            line = raf.readLine();
            if (line == null) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Unxpected End of file while looking for the start of coordinates");
            }
         }

         // --- Reading atoms
         //    1          6      20001001       -0.241185   -0.271881   -0.382124

         while ( (line = raf.readLine()) != null) {
            if (line.contains("-------------------------")) {
               break;
            }

            StringTokenizer st = new StringTokenizer(line, " ");
            if (st.countTokens() != 6) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: coordinate input has 6 tokens\n  Got: " + line);
            }

            GaussianAtom atom = new GaussianAtom();

            // --- Getting element

            // --- Skip 1 token

            st.nextToken();

            int element;
            try {
               element = Integer.parseInt(st.nextToken());
               if (element < 1) {
                  element = 0;
               }
               atom.element = element;
               atom.name = ChemicalElements.getElementSymbol(element);
            }
            catch (Exception ex) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Error while parsing atom: " +
                   (atoms.size() + 1) + " Cannot parse atom number: " + line);
            }

            // --- Getting atomic type

            // --- Skip it for now

            st.nextToken();

            // --- Getting x,y,z

            float xyz;
            try {
               xyz = Float.parseFloat(st.nextToken());
               atom.xyz[0] = xyz;
               xyz = Float.parseFloat(st.nextToken());
               atom.xyz[1] = xyz;
               xyz = Float.parseFloat(st.nextToken());
               atom.xyz[2] = xyz;
            }
            catch (Exception ex) {
               throw new Exception(
                   "parseMolecularGeometryOnly: ERROR: Error while parsing atom: " +
                   (atoms.size() + 1) + " Cannot parse atom's coordinate(s): " + line);
            }

            atoms.add(atom);
         }
      }
      catch (Exception ex) {
         throw new Exception("parseMolecularGeometryOnly: ERROR: " + ex.getMessage());
      }

      return atoms;
   }

   public void parseG03Trajectory(String fileName, int flag_1) throws Exception {
      if (fileName == null || fileName.length() < 1) {
         throw new Exception("Trajectory file name is not set");
      }
      /*
          Time in trajectory (femtosec)    2.000000D-01
          Total energy (au)  -1.175171068D+02
          Total angular momentum (h-bar)   2.904529344D-14
          Coordinates (Bohr)
          I=    1 X=  -4.555420059248D-01 Y=  -5.134698398949D-01 Z=  -7.216904569806D-01
          I=    2 X=  -2.221805396142D+00 Y=  -1.254178689967D+00 Z=  -2.596959873079D-02
          I=    3 X=  -5.952104662883D-01 Y=   1.532677098577D+00 Z=  -8.366401222195D-01
       */

      String line;
      long offset;
      RandomAccessFile raf = null;

      try {
         File f = new File(fileName);
         raf = new RandomAccessFile(f, "r");

         gAtoms = parseMolecularGeometryOnly(raf);
         if (gAtoms.size() < 0) {
            throw new Exception("Didn't find atoms in file");
         }

         /*
           Time in trajectory (femtosec)    2.000000D-01
           Total energy (au)  -1.175171068D+02
           Total angular momentum (h-bar)   2.904529344D-14
           Coordinates (Bohr)
           I=    1 X=  -4.555420059248D-01 Y=  -5.134698398949D-01 Z=  -7.216904569806D-01
           I=    2 X=  -2.221805396142D+00 Y=  -1.254178689967D+00 Z=  -2.596959873079D-02
          */
         while ( (line = raf.readLine()) != null) {
            if (line.contains("Time in trajectory ")) {

               GaussianSnapshot snapshot = new GaussianSnapshot();

               // Add absolute time

               double value;
               String token = "";
               /*
                               try {
                  value = Double.parseDouble(line.substring(30, line.length() - 1).trim());
                  snapshot.addProperty("Time in trajectory (femtosec)", value);
                               }
                               catch (Exception ex) {
                               }
                */

               // parsing            Total energy (au)  -1.175171068D+02
               if ( (line = raf.readLine()) == null) {
                  break;
               }

               try {
                  token = line.substring(18).trim().replaceAll("[D]", "E").replaceAll("[d]", "E");
                  value = Double.parseDouble(token);
                  snapshot.addProperty(TOTAL_ENERGY_KEY, value);
                  propertiesToChart.add(TOTAL_ENERGY_KEY);
               }
               catch (Exception ex) {
                  System.err.println("Cannot parse total energy: " + ex.getMessage() + " got: " + token);
               }

               // Total angular momentum (h-bar)   1.155222129D-13

               if ( (line = raf.readLine()) == null) {
                  break;
               }
               try {
                  token = line.substring(31).trim().replaceAll("[D]", "E").replaceAll("[d]", "E");
                  value = Double.parseDouble(token);
                  snapshot.addProperty("Total angular momentum (h-bar)", value);
                  propertiesToChart.add("Total angular momentum (h-bar)");
               }
               catch (Exception ex) {
                  System.err.println("Cannot parse angular momentum: " + ex.getMessage() + " got: " + token);
               }

               // --- Coordinates (Bohr)
               if ( (line = raf.readLine()) == null) {
                  break;
               }
               double factor = Constants.ONE_BOHR;
               if (!line.contains("Bohr")) {
                  factor = 1.0;
               }

               offset = raf.getFilePointer();
               snapshot.setOffset(offset);

               MolecularGeometry geom = new MolecularGeometry();
               geom.setName(STANDARD_GEOMETRY_KEY);

               // ---Read atomic coordinates
               for (int i = 0; i < gAtoms.size(); i++) {
                  if ( (line = raf.readLine()) == null) {
                     break;
                  }
                  vPoint3f point = new vPoint3f();
                  String buffer = line.replaceAll("[D]", "E").replaceAll("[d]", "E");
                  point.setX( (float) (factor * Double.parseDouble(buffer.substring(11, 32).trim())));
                  point.setY( (float) (factor * Double.parseDouble(buffer.substring(35, 56).trim())));
                  point.setZ( (float) (factor * Double.parseDouble(buffer.substring(59, 80).trim())));
                  geom.addCoordinates(point);
               }

               if (line != null) {
                  snapshot.addGeometry(STANDARD_GEOMETRY_KEY, geom);
                  snapshots.add(snapshot);
               }
            }

            if (line == null) {
               break;
            }

         }

      }
      catch (Exception e) {
      }

      try {
         raf.seek(0);
         raf.close();
      }
      catch (Exception e) {
      }

      logger.info("Found " + snapshots.size() + " snapshots");
   }

   public double[] getAllTerms(String term) {
      if (!propertiesToChart.contains(term) || snapshots.size() < 1) {
         return null;
      }

      List<Integer> references = new ArrayList<Integer> (snapshots.size());
      List<Double> values = new ArrayList<Double> (snapshots.size());
      for (int i = 0; i < snapshots.size(); i++) {
         GaussianSnapshot snapshot = snapshots.get(i);
         try {
            Double value = snapshot.getProperty(term);
            values.add(value);
            references.add(new Integer(i));
         }
         catch (Exception ex) {}
      }
      if (values.size() < 1) {
         return null;
      }
      double[] energies = new double[values.size()];
      for (int i = 0; i < values.size(); i++) {
         energies[i] = values.get(i);
      }
      return energies;
   }

   public int getSnapshotsCount() {
      return snapshots.size();
   }

   @Override
  public float[][] getStructure(int n) {
      return getStructure(n, G03MDTrajectory.TOTAL_ENERGY_KEY);
   }

   @Override
  public float[][] getStructure(int n, String term) {
      if (!propertiesToChart.contains(term) || snapshots.size() < 1) {
         return null;
      }

      if (snapshots == null || snapshots.size() < 1) {
         return null;
      }

      if (debug) {
         logger.info("getStructure: # " + n);
      }

      GaussianSnapshot snapshot = snapshots.get(n);
      MolecularGeometry geom = snapshot.getGeometry(STANDARD_GEOMETRY_KEY);

      if (geom == null || geom.size() < 1) {
         System.err.println(this.getClass().getCanonicalName() + " : no atoms in selected snapshot");
         return null;
      }

      float[][] coords = new float[geom.size()][3];
      for (int i = 0; i < geom.size(); i++) {
         vPoint3f c = geom.getCoordinates(i);

         coords[i][0] = c.getX();
         coords[i][1] = c.getY();
         coords[i][2] = c.getZ();
         if (debug && i == geom.size() - 1) {
            logger.info(i + " : x=" + coords[i][0] + " y=" + coords[i][1] + " z=" + coords[i][2]);
         }

      }
      return coords;
   }

   /**
    * Not implemented yet...
    * @param number int
    * @throws Exception
    */
   @Override
  public void selectStructure(int number) throws Exception {
      throw new Exception(this.getClass().getCanonicalName() + ": selectStructure(int number) is not implemented yet");
   }

   /**
    * Not implemented yet...
    * @param number int
    * @param term String
    * @throws Exception
    */
   @Override
  public void selectStructure(int number, String term) throws Exception {
      throw new Exception(this.getClass().getCanonicalName() + ": selectStructure(int number, String term) is not implemented yet");
   }

}
